Decoherence in quantum spin systems 
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Abstract. Computer simulations of decoherence in quantum spin sys- 
tems require the solution of the time-dependent Schrodinger equation 
for interacting quantum spin systems over extended periods of time. 
We use exact diagonalization, the Chebyshev polynomial technique, four 
Suzuki-formula algorithms, and the short-iterative-Lanczos method to 
solve a simple model for decoherence of a quantum spin system by an 
environment consisting of quantum spins, and compare advantages and 
limitations of different algorithms. 



1 Introduction 

The description of a quantum spin system (below referred to as a central spin 
system (CSS)) interacting with its quantum environment (bath) is among the 
most fundamental problems of theoretical physics. Even if the energy exchange 
between the CSS and the bath is absent (no dissipation) , the system- bath inter- 
action still strongly affects the motion of the CSS due to a loss of phase coherence 
between different eigenstates of the central system. This many-body quantum 
phenomenon is commonly called decoherence. 

Decoherence is fundamental for quantum measurement theory |1I2I3| and for 
condensed matter physics; it can suppress the tunneling of defects in crystals 

spin tunneling in magnetic molecules and nanoparticles 0, or can destroy 
the Kondo effect in a dissipationless manner [Sj. Decoherence is of particular 
relevance for quantum computation since the loss of phase relations between 
different states of the quantum computer may result in an accumulation of errors 
and may prevent the computer from working correctly j7j . A detailed theoretical 
understanding of decoherence would definitely help to alleviate this fundamental 
problem. 

Most theoretical studies of decoherence are based on a model of a single spin 
interacting with a bath of bosons 0] . This model is too simple in the context of 
e.g. quantum computation or tunneling in magnetic molecules. Extensive studies 
of many-spin central systems interacting with other types of environment, such 
as a bath of nuclear spins |S] , are needed. 



A many-spin system interacting with a bath of quantum spins presents a 
fairly complex many-body quantum problem, and numerical simulation is an 
indispensable tool for investigating the long-time dynamics of a decohered CSS. 
One of the most reliable approaches is to model directly the quantum motion of 
the whole system (CSS plus bath) by solving the corresponding time-dependent 
Schrodinger equation (TDSE). For such simulations, the numerical algorithms 
that solve the TDSE should be 1) numerically stable (i.e. conserve the norm 
of the wave function) for all integration times of interest, 2) sufhciently accu- 
rate and allow for controlled increase of the accuracy (e.g. to rule out that the 
loss of phase coherence is due to poor accuracy, rounding errors etc.), 3) effi- 
cient in terms of memory and CPU use, in particular for large spin systems. 
Below we compare three different numerical techniques that have the potential 
to meet these requirements: four Suzuki-formula algorithms 8 9 10 11 12^3], a 
Chebyshev polynomial technique |14I15I16I17| . and the short-iterative-Lanczos 
method 

2 Model and algorithms 

Consider the simple-looking but non-trivial model defined by the Hamiltonian 

L 

H^Jo{Sl + S2f+Y,JnIn -{81+32). (1) 

n=l 

The CSS {81,82), where Si = S2 = 1/2, is coupled to L bath spins {/„} 
{In = 1/2) by a Heisenberg exchange interactions {Jn}- The initial states of the 
spins {In} are assumed to be random and uncorrelated. For the initial state of 
the CSS we take the state with one spin up and the other spin down. We are 
interested in the time evolution of the magnetization of one of the CSS spins, 
e.g. (Sm)- 

A nice feature of the model is that if all Jn — J then, in the large L limit, 
{Sf{t)) can be calculated exactly |21| : 

{Sl{t)) = 1[1 + 2(1 - Ljh^)e-^-''''/^] cos2(Jo - J)t. (2) 

The result Q exhibits an interesting feature: initially, the amplitude of the 
magnetization rapidly decays to zero, then increases again and becomes con- 
stant (1/6) as t — » 00 123. This is similar to the two-step decoherence process 
discovered earlier |22| and can be understood from simple physical arguments 
|21|. The model captures some non-trivial aspects of decoherence, and pro- 
vides a simple test to compare various algorithms for solving the TDSE under 
conditions that are rather demanding from the point of view of algorithmic, 
memory and CPU requirements. We now discuss four different approaches to 
solve the TDSE for models such as 

Exact diagonalization (ED) is the most straightforward approach. Stan- 
dard library routines can be used to compute all eigenvalues and eigenvectors 



of the D X D matrix H {D = 2^+^ denotes the dimension of the Hilbert space 

spanned by the states of the L + 2 spins 1/2). The initial state is represented 

as a superposition of eigenvectors, and the wave function is obtained by 

two matrix-vector muhiphcations of length D and a phase-shift operation on a 

vector. In practice, the amount of memory needed to store the D x D elements 

of the eigenvectors limits the application of this approach to problems with D of 

the order of 10000, which corresponds to systems with about 14 S = 1/2 spins. 

Memory and CPU time of the ED algorithm scale as and respectively. 

Suzuki product-formula algorithms (SP) are based on the approxima- 
tion e-"" « [/2(t) = e-"Hi/2 _ _^-rTHp/2^-tTHp/2 _ ^^-rTHi/2 ^^^^^ ^ 

J2^=i^j- We consider two different decompositions that can be implemented 
efficiently: The original pair-product split-up |8I11| in which Hj contains all 
contributions of a particular pair of spins, and a XYZ decomposition in which 
we break up the Hamiltonian according to the x,y and z components of the 
spin operators |13| . U2(t) is the building block for the fourth-order-in-time ap- 
proximation e~"^ « U4^{t) = [/2(ar)f72(ar)[/2((l — 4a)r)C/2(aT)C/2(aT) where 
a~ 1/(4 — 4^/'^) ^Uj. The error on the wave function is bounded as ||e^'*'^tff(o) — 
U™{t)'F{0)\\ < c„tr" where t = mr and c„ is positive constant. By construction, 
all these algorithms conserve the norm of the wave function and, as a consequence 
are unconditionally stable [5] . These time-stepping algorithms advance the state 
of the quantum system by small time steps r (T[|_ff || <IC 1) and work equally well 
if the Hamiltonian contains couplings to time-dependent external fields ,13;, . For 
a fixed accuracy, memory and CPU time of the n-th order SP algorithm scales 
as D and nt^^^-^^^^ D respectively 

The Chebyshev polynomial algorithm (CP) |14ll5ll6ll7j uses the iden- 
tity<f(i) -limx^oo [Mt\\H\\)I + 2j:^^^Jkit\\H\\)fk{H/\\H\\)\ if'(O). The poly- 
nomials fk{X) are defined by the recursion fk+i{X)<I'{Q) = -2iXfk{X)1'{0) + 
fk-iiX)^{0) for k>l, fo{X)^{0) = If (0), and fi{X)W{0) = ~iX^{0). Using 
standard 14-digit arithmetic, all Bessel functions |Jfe(2)| are zero to machine 
precision ii k > K ^ \z\ + 100 = + 100 and therefore the Chebyshev 

polynomial approximation to \['{t) is accurate to machine precision also (up to 
small rounding errors). Although the CP algorithm is not unconditionally sta- 
ble, it is so accurate that it can safely be used for time stepping (also with very 
large time steps). Note that once t has been fixed, the CP algorithm cannot 
be used to generate reliable information for shorter times. As K is linear in t, 
the computation time required to reach a time t increases linearly with t (and 
D). This linear dependence on t (and the very high accuracy) suggests that the 
Chebyshev polynomial algorithm may be the method of choice if we want the 
solution of the TDSE for a few (very long) times ^Sj- Memory and CPU time 
of the CP algorithm scale as D and tD respectively {K <C D for most problems 
of interest). 

The short iterative Lanczos algorithm (SIL) |15I18I19I2U) is based on 
the approximation e^"^^ « ^-^rPNHPNyjj -^j^ere P/v is the projector on the 
A^-dimensional subspace spanned by the vectors H'F, . . . , H^~^'I'}. We cal- 
culate e"*'^^"^'^™!?' by generating the orthogonal Lanczos vectors in the usual 
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Fig. 1. Left: Magnetization {Sf{t)) as a function of time as obtained by numer- 
ical simulation of two central spins interacting with a bath of L = 10 spins. The 
parameters of model ^ are Jq — 8, Jk — 0.128. Except for the CP algorithm, 
a time step r = 0.05 was used. Right: Comparison of the efficiency of various 
algorithms to solve the TDSE, for the case of the data shown at the left. The 
entry -MP- denotes "machine precision" . CPU times as measured on a Windows 
2000 Athlon XP 1900-F system. 



manner [53], and use exact diagonalization of the resulting N x N tri-diagonal 
matrix for time propagation |15I18I19] . Clearly q-'^'^p^hPn jg ^nj^ary and hence 
the method is unconditionally stable. The accuracy of this algorithm depends 
both on the order N and the state '15"18 19 . In exact arithmetic, e~"^ = 
lim^v^cx) e~*'^^™^'^™!Z', but in practice, the loss of orthogonality during the Lanc- 
zos procedure |22| limits the order TV and the time step r that can be used 
without introducing spurious eigenvalues '23'. Furthermore, we require N <^ D 
because the memory needed to store the eigenvectors (and/or all Lanczos vec- 
tors) is proportional to iV^. In practice, the low-order SIL algorithm may not 
work well if contains contributions from many eigenstates of H with very dif- 
ferent energies, because it is unlikely that all these eigenvalues will be present 
in PnHPn (for small N). Memory and CPU time of the SIL algorithm scale 
as D and N'^Dt/r respectively. In general, N increases with t in a non-trivial, 
problem dependent manner. 



3 Numerical tests 

In Fig. 1, we show a typical simulation result for {Sf{t)), as obtained by the 
CP solution of the TDSE for model The initial fast decay, and subsequent 
reappearance of the oscillations is clearly present. Qualitatively these results 
agree with the analytical (large L) solution (|2Jl. Also shown is the error ||S'jT;]3(i = 
20) — 'P^^t = 20)11 where X is one of the seven algorithms used. It is clear that 
SIL is not competitive for this type of TDSE problem, as already anticipated 
above. The fourth-order pair-approximation is close but still less efficient than 
the CP algorithm, but the other SP algorithms are clearly not competitive. 



The reason that the pair-approximation is performing fairly well in this case 
is related to the form of the Hamiltonian QJ. The present results support our 
earlier finding that the numerical simulation of decoherence in spin systems 
is most efficiently done in a two-step process: the CP algorithm can be used to 
make a big leap in time, followed by the SP algorithm calculation to study the 
time dependence on a more detailed level. From a more general perspective, to 
increase the confidence in numerical simulation results, it is always good to have 
several different algorithms performing the same task. 
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